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Quantum-Classical Transition in Dissipative Double- Well Systems 
— A Numerical Study by a New Monte Carlo Scheme — 

Takeshi MATSUO^^^*, Yuhei NATSUME^ and Takeo KATO^t, 

^ Graduate School of Science and Technology, Chiba University, Inage-ku, Chiba 263-8522 
^Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581 

The thermodynamics of dissipative quantum systems with double-well potentials is studied 
by a path-integral Monte Carlo (PIMC) method without truncation to the two-state model. 
For efficient simulation at low temperatures, we develop a new local update scheme on the 
basis of the approximate decomposition of the Boltzmann weight to Gaussian distributions. 
The localization transition induced by ohmic dissipation is clarified numerically for arbitrary 
potential barriers. 

KEYWORDS: Caldeira-Leggett model, quantum dissipation, path-integral Monte Carlo, localization tran- 
sition 



Quantum tunneling between distinct states subject to 
dissipation has been one of the central problems in both 
physics and chemistry for several decades.^' ^ A number 
of studies on dissipative quantum systems have revealed 
that dissipation strongly affects the quantum nature of 
the system. In particular, the suppression/destruction of 
quantum superpositions has been studied extensively in 
the dissipative two-state model. The most striking fea- 
ture of the two-state system is that localization transi- 
tion occurs at a critical value of dissipation strength for 
ohmic damping. "^"^ This phenomenon can be interpreted 
as a 'quantum-classical transition' driven by the coupling 
with the environment, which is unavoidable for macro- 
scopic objects. 

The dissipative two-state model is useful for the anal- 
ysis of generic systems with a double- well-shaped poten- 
tial in the large-barrier limit. For the small-barrier case, 
however, the truncation to the two-state model cannot 
be justified, and we must analyze the original generic 
system directly. Dissipation effects on generic systems 
are nontrivial, and have been considered only in a few 
studies. Aoki and Horikoshi have discussed localization 
transition in the double-well potential systems by a non- 
perturbative rcnormalization group approach with local 
approximation.^ They concluded that the critical value 
of dissipation strength is larger than that of the two- 
state model. Capriotti et at have studied the thermody- 
namics of dissipative double-well systems by the path- 
integral Monte Carlo (PIMC) method.'' However, they 
have not discussed localization transition, because their 
algorithm based on the conventional update scheme sig- 
nificantly limits the accuracy of the simulation at low 
temperatures. Thus, neither the identification of local- 
ization transition nor the determination of critical dis- 
sipation strength have been achieved so far without ap- 
proximations. 

Recently, a significant improvement in the update 
scheme of the PIMC method has been presented by 
Werner and Troyer for periodic-potential systems.®'^ 
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They introduced a novel update scheme by extending the 
cluster algorithm for classical Ising-spin systems^°"^^ to 
realize efficient change in the paths beyond the poten- 
tial barrier. In addition, they adopted an efficient local 
update scheme by combining the fast random-number 
generation of the Gaussian distributions with stochas- 
tic acceptance. They demonstrated that their algorithm 
is much more efficient than conventional PIMC simula- 
tions,^^' and determined the phase diagram success- 
fully by the accurate calculation of response functions at 
low temperatures. This development of the Monte Carlo 
algorithm, in principle, enables us to study the thermo- 
dynamics of various types of dissipative systems in much 
wider parameter regions. The application of their method 
to double-well systems, however, is not straightforward. 
The cluster update can be implemented successfully to 
double-well systems, whereas their local update scheme 
has difficulty in implementing low acceptance rates at low 
temperatures. This is because their local update scheme 
is based on the fact that a particle in a weak periodic 
potential can move almost freely, whereas a particle in a 
double- well potential is trapped in the whole well-shaped 
structure, even for a small potential barrier. 

In this paper, we propose a new scheme for efficient 
local update in the PIMC simulation, and demonstrate 
its efficiency for large system sizes. By combining this 
improved update with the cluster update by Werner and 
Troyer, we show for the first time the quantum-classical 
transition in generic double- well potential systems with- 
out truncation to the two-state model. We stress that 
our local update scheme has wide application to vari- 
ous classes of models and is not restricted to dissipative 
systems. 

The dissipative system is described by the Caldeira- 
Leggett model:^^'^^ 



The environment properties are uniquely determined by 
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the spectral function 
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In this paper, we restrict our discussion to ohmic dis- 
sipation, for which the spectral function is given as 
J{u}) = m-fuj. The double- well potential is expressed by 
the quadratic-plus-quartic form 

2 



V{x) = 



a 
T 



(3) 



for a barrier height Vq and a spatial distance a between 
two potential minima. In the following discussion, we 
rescalc the length by a/2, the energy by Vo and the fre- 
quency by uiQ = (32Vb/™fl'^)^''^ for a small oscillation 
near the potential minima. Then, in the absence of dis- 
sipation, the system can be characterized by only one 
parameter g — huJo/Vo- The limit g ^ corresponds to 
the semi-classical limit, in which the system can be trun- 
cated to the two-state model using the dilute instanton 
approximation.^ We further define a dimensionless dissi- 
pation strength as a = m-ja^ /{27rfi). By integrating the 
variables of the harmonic oscillators, we can rewrite the 
partition function in terms of path integrals using the 
effective action as 



Z = j Vx{t) exp {-S[x{t)\) , 
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where the memory kernel A;(t) is given as 



k{T) = 



a 



(5) 
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4sin2((7r//35)T)" 

For the PIMC simulation, the path is discretized by the 
Trotter number N, which is taken to be sufficiently large 
to satisfy At = (3/N < 1. After the discretization, the 
action is obtained by the path fragment Xj = x{{(3g/N)j) 
as 

S[{xj}] = -Y,K{i-3)xiXj + AtX^(x2 - If, (7) 

where the long-range kernel including the kinetic energy 
is given by 



AN a 
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for i ^ 0. The value of -fl'(O) is taken so that the 
zero-frequency Fourier component of the kernel becomes 
zero. In the conventional update scheme, one of the 
path fragments is changed as x^"^ = xf'^ -|- e for a 
fixed j. Then, this change is accepted with the prob- 
ability min(exp(-(S'[{xf ^}] - S[{xf'^}])), 1). Here, the 
displacement e is generated randomly by the uniform dis- 
tribution with the range [—5, S], where S is chosen appro- 
priately to obtain the optimal acceptance rate. This al- 
gorithm, however, takes considerable time for one sweep 
of local updates on all the path fragments. 
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Fig. 1. Boltzmann weight as a function of position variable x 
for double-well potential (solid line) and approximated poten- 
tial (broken line), which is a sum of four Gaussian distributions 
(other lines). 

We now introduce a new update scheme. First, the 
Boltzmann weight of the potential term is approximately 
decomposed into M pieces of the Gaussian distribution 
as 



M 



g-A.y(x)^^^^^-A(.- 



(9) 



where A, (Jm and arc determined by fitting. For the 
double- well potentials, V{x) — (1 — x^)^, and this ap- 
proximate decomposition can be carried out with a small 
value of M. One example of fitting is shown in Fig. 1 
for At = 100/256 ~ 0.4. Here, we used M = 4 Gaus- 
sian distributions, and the parameters were obtained as 
A = 3.7683, -(Ti = (74 = 0.9068, -a2 = 0-3 = 0.5858, 
^1 = /i4 = 1.1169 and /i2 — ^J.3 — 0.3963. As seen 
in Fig. 1, the approximate fitting is almost in agree- 
ment with the original weight. In our calculation, we 
checked that the deviation from the exact weight is not 
distinguishably different for the obtained results up to 
TV = 1024. By this decomposition, the partition function 
is written in the form 



^ = E / ll<ixjW{{x^},{r 

{rrij} j 



(10) 



Here, the new variables {mj} play the role of auxiliary 
fields. Our update scheme consists of two steps. In the 
first step, the variables {nij} are chosen for fixed {xj} 
with the probability t{mj \xj). For this update, the weight 

is modified as VF({a:j}, {rrij}) = f{{xj}) x Yl^ w{xj,mj), 

where w{xj,mj) = exp(-A(a:;j -cTm,)^) and f{{xj}) 
is independent of {rrij}. Then, the transition probability 
is given by the extended detailed valance condition^^ as 

w{xj,mj) 



t{mj\xj) = 



Em, w{xj,mj)' 



(11) 



Note that this update needs only the local informa- 
tion of each path fragment, thus, can be performed ef- 
ficiently. In the second step, the path fragments {xj} 
are updated for fixed {mj}. In this update, it is cru- 
cial that for fixed {mj}, the effective action defined by 
^{{xj}j {n^j}) = exp(— 5) is written in a quadratic form 
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of {xj}. Using a Fourier transformation, the effective ac- 
tion can be written in the normal form 

^ = ^(c„|5;„|2+AAf|x„-a„|2) - ^In/im., (12) 

n j 



, nn n a , , 
c„ = -^sm-- + — H, 



(13) ^ 



where Xn and (T„ are the Fourier components of Xj and 
CTmj, respectively. The effective action can be modified 
further as 

2 



S = l](Cn + XN) 



\Nd„ 



Cn + XN 



+ g{{mj}), (14) 



where giirrij}) is independent of {xj}. Then, the new 
path is efficiently generated by the Gaussian distribu- 
tion with mean \N(7n/{cn + XN) and standard devia- 
tion (c„ -I- AA^)^^/^/2 for both real and imaginary parts 
of .T„.^^ Note that all the path fragments are updated 
immediately without rejection. 

In the presence of large potential barriers, global 
change of the path across the barrier is strongly sup- 
pressed, even in the improved local update. To overcome 
this difficulty, we adopt the cluster update by Werner 
and Troyer.'^ We first focus on the mirror symmetry 
V{x) = V{—x) for the potential energy, and flip the sign 
of the path without varying its amplitude, \xj\. Since the 
potential term in eq. (7) is unchanged by this update of 
the path, we can disregard this term in the following dis- 
cussion. We define the Ising variables as Xj = Sj\xj\, and 
express the effective action by these variables as 



'^'[{Sj}] — — ^ JijSiSj, 



(15) 



where Jij 



K{i — j)\xi\-\xj\. This is simply the one- 
dimensional long-range Ising model. Then, the Wolff 
cluster algorithms can be applied directly. We note that 
the cluster update can deal with quantum tunneling pro- 
cesses efficiently through the; creation (annihilation) of 
instantons and the change in their positions. 

We show the effectiveness of our improved local up- 
date in Fig. 2. Here, the CPU time in a improved cal- 
culation is measured in a cluster computer consisting of 
eight Pentium 4 processors for the Monte Carlo sampling 
within the autocorrelation time in the correlation func- 
tion {x{0)x{P/2)). The main figure in Fig. 2 shows the 
comparison between the conventional algorithm and the 
new algorithm for the local update at a = (no dissipa- 
tion) without the cluster update. As seen in this figure, 
our algorithm for the local update shows significant im- 
provement in the simulation time. The improved CPU 
time is approximately proportional to N for large sys- 
tem sizes, for which the autocorrelation time is almost 
independent of iV. This system-size dependence indicates 
that the efficiency of our simulation is limited by the 
Fourier transformation, which requires time on the order 
of NlogN. The inset in Fig. 2 shows the CPU time mea- 
sured for a > (nonzero dissipation) using both the new 
local update and the cluster update. This result shows 
that the efficiency of our algorithm is less affected by 
dissipation strength a. 
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Fig. 2. Effective simulation time for Monte Carlo sampling 
during autocorrelation time measured in correlation function 
(a;(0)a;(;8/2)). Main panel: comparison between conventional and 
new local update scheme with no cluster update for a = (no 
dissipation) and g = 2. Inset: a-dependence of effective simula- 
tion time for both cluster update and new local update. 




Fig. 3. Calculated susceptibility of double-well potential system 
as a function of damping parameter a. The curves correspond to 
the Trotter numbers N = 2^^, 2^^, and 2^^ from the bottom 
to the top. The inset shows the scaling fitting for N = 2^"^ ~ 
1.3 X 10^ in the region 1.42 < a < 1.62. 

Next, let us move onto the localization transition 
in double- well potential systems. We expect that this 
transition belongs to the Kosterlitz-Thouless type be- 
cause of the same symmetry as the long-range Ising 
model with inverse square-law interactions.^^ In accor- 
dance with ref. 19, the transition is determined by 
calculating the uniform susceptibility defined as x = 
N(^^^jXiXj/N'^)/{Y^^x1/N). We use renormalization 
group analysis^'^ suggesting that 

X ~ eM^ia-a,)-'' + B{a-a,y + 0{{a-a,y)]. (16) 

The critical value Uc is then obtained by fitting. We show 
the data of the susceptibility at (/ = 2.0 for large system 
sizes up to A?' = 2^^ ~ 1.3 x 10^ in the main panel of 
Fig. 3. To plot one point in this figure, we used 64 in- 
dependent runs, each of which consisted of 500 Monte 
Carlo steps (MCS) for thermalization and next 600 MCS 
for measurement. In all the runs, we checked that the au- 
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Fig. 4. Phase diagram of dissipative double-well systems in p-a 
plane, where a is the dissipation strength and g = huio/Vo is 
the zero-point quantum fluctuation normalized by the potential 
barrier. The broken line is drawn to guide the eyes. 

tocorrelation time measured in the correlation function 
(x(0)x(/3/2)) is much less than 300 MCS. The inset shows 
the fitting for 1.42 < a < 1.62, where the data has almost 
converged on one curve with respect to the change in the 
system size. Stable fitting is obtained for each value of 
g, and the critical dissipation strength ac is successfully 
determined within reasonable statistical errors. 

Finally, the phase diagram is given in Fig. 4, where 
the closed squares show the critical values of dissipation 
strength for localization, We stress that it is the first 
time that transition points in the present system have 
been determined as a function of g without approxima- 
tion.^^ The error bar shown for the data a.t g = 1.5 ex- 
presses the statistical error estimated from fitting using 
five independent Monte Carlo results. The error is re- 
duced for larger values of g, and becomes smaller than 
the symbol size at g = 4. As g increases, is enhanced 
from the value ac = 1, predicted in the two-state limit 
5f — > (shown in the figure by the open square). Further- 
more, the tendency of saturation is clearly observed for 
larger values of g. The behavior of ac indicates that even 
for the infinitesimal potential barrier (Vq <C tuvo), a finite 
strength of dissipation is sufficient for the localization 
transition. This impressive result has already obtained 
for dissipative periodic-potential systems, where the crit- 
ical value of dissipation strength is obtained as ac = 1 
in the small-barrier limit Vq/Hujo 0.^'^^ Note that the 
periodic-potential system is special because of the exis- 
tence of a property called 'duality', which strongly sug- 
gests that the critical value ac is independent of the po- 
tential barrier. Actually, for the large-barrier limit, the 
critical value is also given as ac = 1. Because double- 
well potential systems do not have such a property, the 
critical value can be derived from the limiting value. 

The improved local update developed by us can be ap- 
plied to other models, one of which is the nonlinear har- 
monic oscillator system such as the chain. The present 
scheme is also useful for analyzing the nonlinear effects 
of phonons in dielectric substances. The extension of our 
algorithm to these systems will be studied elsewhere. 
The idea of Gaussian decomposition also suggests that a 



model with continuous variables can be mapped approx- 
imately into another model only with discrete variables 
after integrating out the continuous variables. This pro- 
cess will give a new insight into the original continuum 
model. 

In summary, we considered the localization transition 
in dissipative double-well systems by a PIMC simula- 
tion improved by an approximate decomposition of the 
weight by several Gaussian distributions. We determined 
the transition point by fitting to the scaling expression 
for the Kosterlitz-Thouless transition. For potential bar- 
riers that are small compared with quantum fluctuations, 
the critical value of dissipation strength becomes larger 
than 1, and saturates at a finite value. We expect that the 
decomposition into Gaussian distributions is helpful for 
both creating algorithms and analyzing of other models 
with continuous variables. 
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